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Abstract — THIS PAPER IS ELIGIBLE FOR THE STUDENT 
PAPER AWARD. 

This paper extends some approximation methods that are 
used to identify closed form Bit Error Rate (BER) expressions 
which are frequently utilized in investigation and comparison 
of performance for wireless communication systems in the 
literature. By using this group of approximation methods, some 
expectation integrals, which are complicated to analyze and have 
high computational complexity to evaluate through Monte Carlo 
simulations, are computed. For these integrals, by using the 
sampling property of the integrand functions of one or more 
arguments, reliable BER expressions revealing the diversity and 
coding gains are derived. Although the methods we present are 
valid for a larger class of integration problems, in this work we 
show the step by step derivation of the BER expressions for a 
canonical cooperative communication scenario in addition to a 
network coded system starting from basic building blocks. The 
derived expressions agree with the simulation results for a very 
wide range of signal-to-noise ratio (SNR) values. 

I. Introduction 

There has been a great effort for mitigating the detrimental 
effect of fading in wireless communication systems for a long 
time. Over the past few years, cooperative communication 
techniques, which rely on the spatial diversity thanks to 
the overhearing nodes within the network, have been in the 
center of these efforts [1|. In achieving this cooperation, 
various methods making use of relay nodes have been pro- 
posed (21, 0. The simplest method with demodulate- and- 
forward (DMF) type relays has an inherent error propagation 
problem which must be handled at the destination node. For- 
tunately, when the destination uses the maximum-likelihood 
(ML) detection [ 3 ] or cooperative maximal ratio combining 
(C-MRC) |4] technique, this problem is overcome and the 
maximum diversity order is achieved in the system. The 
canonical cooperative communication scenario, where a single 
relay node assists the transmission from the source node to the 
destination, is investigated for the destination with the C-MRC 
operation in E). The main advantages of C-MRC include its 
reduced computational complexity with respect to the optimal 
ML method for increasing number of relays with practically 
no loss of BER performance and simpler analysis. 



For the end-to-end BER analysis of C-MRC, in |4], the 
diversity order is considered, which is obtained following a 
series of detailed upper bounding operations. The closed form 
BER expressions including the coding gain terms in addition 
to the diversity order terms would be of importance for this 
canonical scenario and the one with multiple relays. 

Recently in (5), (6J, an approximation method applicable to 
expectation integrals of a single random variable is proposed 
for the analysis of fading channel BER performance of coop- 
erative systems. The crucial observation of these papers is that 
the Q-function appearing in the instantaneous BER function 
can be well approximated with a Dirac delta well especially 
in the high SNR region. However, the upper bound utilized on 
Q-function to obtain the impulse location leads to degradation 
in the approximation results. More importantly, the analysis 
lacks rigorous handling of the properties of the sampling 
function at points other than the critical point. Moreover, 
the analysis in 0, 01 takes directly the Q-function as its 
sampling function without any consideration on the SNR value 
for which the approximation is made. In [7 |, the SNR regions 
are identified using the relative error of each possible selection 
of the sampling function, however a sound analysis of the 
sampling property is still not provided. 

In this work, we have four goals: (i) A more rigorous 
method to show the applicability of the sampling property, 
(ii) Characterization of the SNR regions for sampling property 
which yields improved approximate expressions, (iii) A more 
accurate identification of the impulse function parameters, (iv) 
Generalization of the idea of sampling to integrals involving 
more than one variables (here for the canonical cooperative 
scenario, 2 variables) to obtain closed form BER expressions 
approximating these integrals. Finally, as an application of 
sampling property, we sketch the method to obtain closed form 
expressions that well- approximate the simulated performance 
of a sample network coded scenario from [[8]. 

II. Canonical Cooperative Communication Model 

and Instantaneous End-to-End BER Expression 

The canonical relay network includes the link from the 
source node S to the relay node R called S — R, the one from 



R to destination D named R — D in addition to the direct 
link from S to D shown S — D [4]. The complex channel 
fading coefficients related to these links are shown in Fig. Q] 
In the first time slot only S transmits the data symbol x and 
in the second slot only R transmits its own possibly erroneous 
detection result x R . For the sake of simplicity in analysis, as 



the instantaneous end-to-end BER expression is found as 0] 




RD 



Fig. 1. Canonical cooperative communication model 



in H, we assume x e {+1,-1}. In both of the transmissions 
the average transmit power is taken as P. The independent 
channel fading coefficients hsR, h R o, and Hsd are assumed 
to follow zero mean circularly symmetric complex Gaussian 
probability distributions such that hij ~ CiV(0, of ■), where 
(ij) e {SR, RD, SD}. Then the received signals at R and D 
following the first slot are written respectively as 

VSR = hsRX + nsR and ysD = h S DX + n S D, (1) 

where usr and usd denote the independent white complex 
Gaussian noise terms at R and D with identical distribution 

CN(0,N ). We define 7 = ^, fU 



then the instantaneous SNR 
values for S-R and S-D links are jsr = j\hsR\ 2 and 
Isd = tI^sdI 2 - These instantaneous values are exponentially 
distributed with respective expectations 7ct|d and *fcr 2 SD . In 
accordance with the DMF method, prior to the transmission 
in the second slot R should detect x using the optimal (ML) 
rule as follows. 



x R = argmin \y SR - h S Rx\ , 
xG{ + l,-l} 



(2) 



where we assume perfect knowledge of hsR at R. 

The observation of D after the second time slot is then 
Vrd = h RD x R + n RD , where n RD ~ CN(0,N ) and the 
corresponding instantaneous SNR j R d = j\Ii R d\ 2 is also ex- 
ponentially distributed. The aim of D is to correctly detect the 
data symbol x of node S by making use of all instantaneous 
channel gains. If we follow the C-MRC technique to linearly 
combine two observations at D as suggested in [4], we obtain 
the following detection result. 

x D = argmin \wiy S D + w 2 yRD ~ (wih S D + w 2 h RD )x\ 2 , (3) 
xG{+l-l} 

where w\ is the weight coefficient corresponding to the S — D 
link and selected to be equal to hg D just as in the well- 
known MRC case with no relaying operation. On the other 
hand, for the observation urd, which corresponds to the 
relayed communication over the links S — R and R — D, 
the coefficient w 2 should be re-defined to reflect the possible 
error propagation on these two hops. Authors, in [4 |, propose 
a single equivalent channel for representing these hops and 
assign the equivalent instantaneous SNR of this channel as 
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(4) 



±{Q-([l 
where Q(x) = 

and P h RD = Q(v / 2t^). If we further define w 2 = ^h* RD , 



[i-Prd]Psr)} /2 
exp(-z 2 /2)dz, P b SR = Q(V2J^ 



■Q(y/2^R))Q 



V2(lSD +7eq) 
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(5) 



Clearly, it is analytically hard to evaluate a triple integral for 
the instantaneous BER given in (0) over related distributions 
of the random variables jsr, Isd, and jrd (and j eq which is 
a complicated function of jsr and jrd) to reach the average 
BER value. In [4], only the diversity order for the BER- 
SNR curve is identified following a set of upper bounding 
techniques and the result is 2 since two transmissions are 
made on independent paths in the network and D accounts 
for possible relaying errors. In Section IIII-EI we present novel 
closed form expressions for the average BER of this system 
such that the coding and the diversity gains are identified 
separately. 

III. Sampling Property of the Q-Function for 
Generalized Expressions 

This section starts with the basic result of [6| with a slight 
improvement on obtaining the critical point of the related 
sampling function. We also exemplify the problem with this 
technique for the low SNR region and propose a method to 
identify threshold value of SNR over which the technique 
is applicable. We further generalize the sampling property to 
obtain good approximations for lower SNR values. Finally, we 
generalize this sampling property to integrals involving more 
than one variables. 

A. Basic Problem and its Solution 

Assume that the following expectation integral of an instan- 
taneous probability of error function Q(y/X) is to be evaluated 
for a random variable X with pdf fx{x) 0- 



h = E x {q(Vx)} = J™ Q(Vx~)h 



(x)dx. 



(6) 



After the change of variables operation x —> t N , we have 
the integrand Q(VW)Nt N ~ 1 fx(t N ). Here we define the 
following constituent functions of the integrand. 



q(t) 4 Q(Vt N ), c(t) 4 Nt N -\ f{t) 4 f x (t N ). (7) 

In 0, h(t N ) = q(t)c(t) is defined and claimed to be a 
unimodal function of t with a critical point satisfying = 2. 
However, this analysis does not show that h(t N ) takes value 
of zero for any other value and on the contrary, when any 
other finite t N value is inserted in Eqn. (52) of it is easy 
to show that h(t N ) assumes infinity. Moreover, the critical 
point = 2 is obtained after applying an upper bound on 
the Q-function. Let us start with the analysis of the integrand 
I(t) = q(t)c(t)f(t) for three distinct regions of t. 



I = i(t)dt + / I(t)dt + / I(t)dt. 
Jo Ji- Ji+ 



(8) 



It should be emphasized that the value Io in ® is independent 
of N and hence we may investigate the behaviour of I(t) 



asymptotically (for N — >> oo). Moreover, in this work, we take 
f(t) = s^tr^satr (Rayleigh fading assumption) although 



Combining CE2> and fl3 for t > 1 and SiVi? > 2 (in dB 
scale roughly for values larger than 3 dB) we get 



the following steps can be generalized to other pdfs. Initially, 
for the < t < 1 region, we have 
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Using O we reach the result that for SNR > 2, h(t N ) 
(including the Q-function) is better represented by a Dirac 
delta with respect to g(t N ) (including the exponential pdf). 
Also, for x > 1, one can show that 



0, 



(9) 



where the L'Hopital rule is applied for liniAr^oo c(t) term. As 
a result, the first integral in © evaluates to asymptotically. 
Next, for the t > 1 region, we have 



: Q(y/X) > 



T ex P 

3 



exp (-< 



SNR 



fit) 



(15) 



lim lit) 



= ( lim q{tj) ( lim c(t)/(t)) 

= (0) 



taking S^Vi? < 1/3 (roughly less than —5 dB). Consequently, 
for lower SNR values, g(t N ) fits better to sampling function 
definition. Firstly, the position of the Dirac delta that approx- 
imates g(t N ) can be obtained by finding the critical point t£ . 
We equate the first derivative of g(t N ) with respect to t to 0: 



lim 

iV— >-oo 
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whose solution is 



SNR 1 eX P 



\SNR J 



N — 1 
N 



= 0, 



SNR. 
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Eqn. (TTrb gives us the asymptotic critical point ^ 



(16) 



(17) 
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(10) 



following application of the L'Hopital rule twice. Hence, in the 
asymptotic sense, the third integral in © does not contribute 
to the result as well. Therefore, we obtain 



Io = f lim I(t)dt, 

A- 



' t N exp 

(ms) J limyv^oo t% = SNR. Secondly, the weight of the correspond- 

ing Dirac delta is found as 1 due to the normalization property 
of the pdf. Hence for SNR < 1/3, we may use the approxi- 
mation I w J °° Q(Vx~)5(x - SNR)dx = Q (y/SNR^. 

For SNR > 2, we write I « / °° cS{x - t?)f x (x)dx = 
exp (^—^^j, where the impulse weight is found using 



(11) 



which shows us that the integrand I(t) may be well- 
approximated by a Dirac delta at t = 1 for N — » oo. Moreover, 
it can be shown that this Dirac delta approximation also holds 
for the functions h(t N ) = q(t)c(t) and g(t N ) = f(t)c(t). 
Herein we are going to assess the limiting SNR value above 
which h(t N ) can be safely approximated by a Dirac delta. 
In order to obtain this threshold, we are going to compare 
convergence rates of h(t N ) and g(t N ) in Section Ull-Bl 

B. Rates of Convergence for Constituent Functions 

If we remember d8}, two regions that asymptotically yield 
no contribution to the integral should be examined for conver- 
gence rate comparison. For the < t < 1 region, both q(t) 
and f(t) converge to nonzero constants in the limit. On the 
other hand, the function c(t) = N^ -1 converges to 0, which 
means that for large N, h(t N ) and g(t N ) converge to with 
the same rate. Hence we will concentrate on the other region: 
t > 1. In this region we are going to compare h(t N ) and g(t N ) 
starting with the Chernoff upper bound on the Q-function. 



SNR 

the alternative definition of Q-function as c = / °° h(t N )dt 
f °° Q(^/x)dx = \. We propose a simple alternative to 
the method in [6] to evaluate the location of the impulse 
approximating h(t N ). For a sufficiently large value of N, we 
pose finding the critical point as an unconstrained optimization 
problem and employ numerical search to find the solution. In a 
narrow neighbourhood of t = 1, we search for the critical point 
of h(t N ) and we find = 1.4157 just after 4 iteration steps, 
whereas in =2 was considered. The approximation 
methods as well as the simulation result are given in Fig. |2] 
According to Fig. O the methods approximating h(t N ) as a 



q(t) 



= Q(V5)<|exp(-|). (12) 



In addition, for t > 1 and SNR > 2 it is easily shown that 
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_Q-Apprx. £'= 1.4157 
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Fig. 2. Approximating the integral Iq using various methods 



Dirac delta (one proposed in [6| with square markers, and the 
one we propose with circular markers) are quite consistent for 
high SNR values. However the method using unconstrained 
optimization for searching impulse location is the better one 
with close approximation for SNR > 3dB as detailed in 
equation (fT4l) . For low SNR values, on the other hand, only the 
method selecting g(t N ) as the sampling function (plus shaped 
markers) is in a sensible proximity to the simulation result. 
This shows us that for a close approximation of the integral 
Io over the whole SNR region, we need to use a piecewise 
function as a result of (ITU) and (TT5T) . In the following sections, 
with more than one variables case, we are going to use only 
the sampling property for the Q-function since we are mostly 
interested in high SNR performance. 

C. Two-Variable Sampling Property 

In this section we are going to base our discussion on the 
following integral involving two variables in the integrand. 



D. Sampling in Single Dimension for Functions of Two Vari- 
ables 

Unfortunately, not every integrand function can be simply 
appr oxima ted with a two dimensional Dirac delta as in Sec- 
tion IIII-CI As an example, the instantaneous BER function 

min{x, y}) can be shown to diverge for some points 
other than (t,u) = (1,1) following the suitable change of 
variables operation. However, it is still possible to analyze this 
function using the fact that Q{y2 min{x, y}) < Q(y/2x) + 
Q(y/2y) and the sam pling p roperty for single variable func- 
tions given in Section llTl-AI Then we can reach the following 
approximation for this expectation integral, which is shown to 
perfectly fit the simulation result in Fig. [4] 
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Q^axx + a 2 y)fx(x)f Y (y)dxdy, (18) w i m p 1 = (i 



E. BER analysis for the Canonical Cooperative Model 

The end-to-end instantaneous BER function in ([5]) can be 
written as the sum of two terms: P b = Pi + P 2 . Let us start 

V2{lsn+l eq ) 1 whichisa 



where a\ and a 2 are positive constants. Similar to 
the single dimension analysis we define h(t N ,u N ) = 

Q( + a 2 u N )N 2 t N - 1 u N - 1 based on the Q-function. 
Here, it is easy to generalize the asymptotic analysis for 
h(t N , u N ) with N —> oo to show that it is well- approximated 
by a two-dimensional Dirac delta at (t,u) = (1,1). This is 
further exemplified in Fig. [3] for N = 1000 and a\ = a 2 = 2. 
Such a higher dimensional generalization for sampling prop- 
erty is not made in 0-GI, but whenever a two-dimensional 
integrand is encountered, the double integral is approximated 
by two single variable integrals yielding coding gain errors 
in the final expressions. Through unconstrained optimization 



function of three v ariable s, "ysR, Jrd, and jsd- Similar to the 
analysis in Section Hll-D[ function Pi can not be approximated 
with a Dirac delta directly, due to the variable j eq defined over 
the instantaneous SNR values of the two-hop link, 7^ and 
jsd- Therefore, we define two terms that are asymptotic in 
Ird and jsd following the approach in Eqn. (42) of 



plRD 



= lim Pi 



■Q(V*ysr))Q 



V2(jSD + IS r) 



(21) 




P^flA lim P 1 =Q \y^ys~D~+7RD)} 

to approximate Pi with the sum of these two terms. In this 
way, Pi is now the sum of two functions both of which 
have two arguments and are suitable for an approximation 
with impulse functions. It should be noted that in the ap- 
proach utilized in Section IIII-D1 Q(y/2x) and Q(\/2 y) are 
also asymptotic terms. Using the result of Section IIII-CI 
approximate expectation of Pi is evaluated to be 

poo POO 

h~ / Pl RD f 1S R(lSR)f 1SD (lSD)&lSR&lSD 

Jo Jo 
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(22) 



Fig. 3. Function h(t N , u N ) for N = 1000 

solution, the critical point of h(t N ,u N ) is computed as 
(tf , ) = (0.8197, 0.8197) and the weig ht of the Di rac delta 
is analytically found as c = J °° J °° Q(y/aix + a 2 y)dxdy = 
4q ^ Q2 . As an example for a\ = a 2 = 2, the approximation for 
Ii is 



Defining similar asymptotic terms for P 2 , we reach 

P OC POO 

J 4 « / / P2 RD flSR < C1SR)f 1S D < C1SD)&lSR&lSD 

Jo Jo 



■ exp 



1.7564 + 1.3737 \ 



(23) 
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■ exp 



2(0.8197) \ 
SNR ) ' 



(19) 



16SNR 2 

The result in (|T9l) is in accordance with the average BER 
analysis result of the MRC technique applied on two parallel 
branches and is also very close to the simulation result for 
mid- to high SNR values as given in Fig. [H 



4SNR 2 V SNR 

Finally, summing the results of (l22t and (l23t we obtain and 
plot the approximate expectation of P b as Is + 1 \ in Fig. [4] 
together with the simulation result. It is seen that the analysis 
proposed in this work yields an extremely good approximation 
to the end-to-end average BER of the canonical cooperative 
communication system by giving the closed form expression 
as a product of the coding and diversity gain terms. 
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Fig. 5. Approximating the average BER for the network coded system 



Fig. 4. Approximating integrals I±, I2, and 1 3 + I 4 



IV. BER Analysis for a Network Coded System 



In this section, we are going to sketch the BER analysis 
for a sample network coded cooperative system (given by 
Eqn. (1) and Fig.l of (8|) with single destination node and 3 
source nodes that transmit network coded data in 4 orthogonal 
transmission slots. This DMF-type relaying system model is 
essentially the same with the canonical cooperative model 
of (H with the additional XOR operations at the intermediate 
(relay) nodes and extra time slots. Utilizing the equivalent 
channel approach for the network coded data, one may reach 
the following approximate instantaneous probability of error 
function for the data bit of node l,i.e., u\. 



P(ui / Ul) « (1 ~Pe4)Q 
( 



+Pe4Q 



V2 (\h!\ 2 - \h eq4 \ 



{ V |/il|2 + 7^ J 



( \ 

V^(|fel| 2 + |/leg4| 2 ) 
\ 

+ Q^2(|/ii|2 + |^ 2 |2) 



(24) 



where p e 4 is the probability of error for the transmitted 
network coded data in the fourth time slot, |/ii| 2 , |^2| 2 , and 
I J14 1 2 denote the independent instantaneous SNR values for the 
transmissions in the corresponding slots from the intermediate 
nodes to the destination node and follow exponential distribu- 
tions with mean value equal to SNR. Similar to the equivalent 
channel definition for the canonical cooperation model, | h eq 4 | 2 
is defined to represent the two-hop transmission that ends in 
the fourth time slot. We imm ediatel y make use of the sampling 
properties derived in Section Ull-Cl to approximate the integral 
in (1241) and after a few steps we obtain 



h 3 exp r_ 2 ( Q - 8197 ^ 

8SNR 2 P V SNR J 



16SNR 2 eXP 



1.3049 \ 
" SNR ) 
3.1301 



SNR 



(25) 



One can distinguish the diversity order of 2 and the coding 
gain using the closed form approximation (l25t . Following the 
same procedure for other two nodes' data bits, we plot Fig. [5] 
in order to visualize the approximation results for the network 
coded system. Even for low SNR values, the proposed Dirac 
delta approximation for multi-dimensional integrals is in a 
good agreement with the simulation results. 



V. Conclusion 

In this work we rigorously investigate the applicability 
of one dimensional sampling property to some frequently 
encountered expectation integrals using asymptotic analysis 
for distinct regions of the argument of the integrand functions. 
Identification of the threshold SNR value for selection of the 
most suitable part of this integrand gives us the opportunity 
to define piecewise approximations that are very close to 
the simulation results over the whole SNR region. More- 
over, we extend our sampling property analysis to higher 
dimensions and through some basic examples, we build up 
an approximation method to give closed form expressions 
for the performance evaluation of a canonical cooperative 
communication scenario. These extensions are further shown 
to cover a sample network coded scenario yielding closed form 
expressions that are very close to the simulation results. 
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